Add derivatives for `besselix`, `besseljx`, and `besselyx` #350

devmotion · 2021-09-28T23:08:32Z

This PR adds derivatives for the non-holomorphic functions besselix, besseljx, and besselyx.

codecov · 2021-09-29T01:19:13Z

Codecov Report

Merging #350 (a9f9f52) into master (0c41812) will increase coverage by 0.19%.
The diff coverage is 100.00%.

@@            Coverage Diff             @@
##           master     #350      +/-   ##
==========================================
+ Coverage   92.75%   92.94%   +0.19%     
==========================================
  Files          12       12              
  Lines        2706     2765      +59     
==========================================
+ Hits         2510     2570      +60     
+ Misses        196      195       -1

Flag	Coverage Δ
unittests	`92.94% <100.00%> (+0.19%)`	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

Impacted Files	Coverage Δ
src/chainrules.jl	`100.00% <100.00%> (ø)`
src/beta_inc.jl	`90.37% <0.00%> (-0.31%)`	⬇️
src/gamma.jl	`94.90% <0.00%> (+0.08%)`	⬆️
src/erf.jl	`100.00% <0.00%> (+2.04%)`	⬆️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 0c41812...a9f9f52. Read the comment docs.

devmotion · 2021-10-01T08:03:59Z

Since it's the first time I implement complex derivatives of non-holomorphic functions it might be good if someone more experienced could have a look at the PR. Is there anything that could be improved @oxinabox? It seems the real(conj(..) * ...) could be optimized in the same way as ChainRules._realconjtimes. I was wondering if one could move https://github.com/JuliaDiff/ChainRules.jl/blob/master/src/rulesets/Base/utils.jl to ChainRulesCore and make realconjtimes and imagconjtimes available as part of the official API.

oxinabox · 2021-10-01T10:29:01Z

Since it's the first time I implement complex derivatives of non-holomorphic functions it might be good if someone more experienced could have a look at the PR.

@sethaxen perhaps?

sethaxen · 2021-10-01T10:42:54Z

Great! I'll review tonight or tomorrow.

devmotion · 2021-10-01T12:48:08Z

I was wondering if one could move https://github.com/JuliaDiff/ChainRules.jl/blob/master/src/rulesets/Base/utils.jl to ChainRulesCore and make realconjtimes and imagconjtimes available as part of the official API.

Do you think I should open an issue or PR in ChainRulesCore regarding this question?

sethaxen · 2021-10-01T13:15:22Z

I was wondering if one could move https://github.com/JuliaDiff/ChainRules.jl/blob/master/src/rulesets/Base/utils.jl to ChainRulesCore and make realconjtimes and imagconjtimes available as part of the official API.

Do you think I should open an issue or PR in ChainRulesCore regarding this question?

Yes, I think that's a good idea.

devmotion · 2021-10-01T15:16:34Z

Yes, I think that's a good idea.

I made a PR: JuliaDiff/ChainRulesCore.jl#474

sethaxen

Overall looks good. Most of my comments on besselix apply also to besselyx and besseljx.

sethaxen · 2021-10-01T21:49:34Z

src/chainrules.jl

+    project_x = ChainRulesCore.ProjectTo(x)
+    function besselix_pullback(ΔΩ)
+        ν̄ = ChainRulesCore.@not_implemented(BESSEL_ORDER_INFO)
+        a = (besselix(ν - 1, x) + besselix(ν + 1, x)) / 2


One can use the recurrence relations to write this as (besselix(ν ± 1, x) ± besselix(ν, x)) * ν / x, which allows to reuse Ω. It would just require some special-casing to handle x=0 gracefully.

I just reused the derivatives that are used for besseli etc. They were defined in this way in ChainRules originally. When I copied them to SpecialFunctions I wondered about the motivation for choosing https://functions.wolfram.com/Bessel-TypeFunctions/BesselI/20/01/02/0003/ over https://functions.wolfram.com/Bessel-TypeFunctions/BesselI/20/01/02/0001/ or https://functions.wolfram.com/Bessel-TypeFunctions/BesselI/20/01/02/0002/. I assumed the currently used definition is simpler since it does not require to handle x = 0 in a special way.

I think it might be best to use the same relations for besselix etc. as for besseli etc. and, if desired, change them to a different form in a separate PR.

That's right, and they're the same in DiffRules as well. If you like, you can keep them similar to besseli, etc for now, and a future PR could update all of the rules.

sethaxen · 2021-10-01T22:08:34Z

src/chainrules.jl

+        ν̄ = ChainRulesCore.@not_implemented(BESSEL_ORDER_INFO)
+        a = (besselix(ν - 1, x) + besselix(ν + 1, x)) / 2
+        b = - sign(real(x)) * Ω
+        x̄ = project_x(conj(a) * ΔΩ + real(conj(b) * ΔΩ))


It's more efficient to compute -sign(real(x)) * real(conj(Ω) * ΔΩ)) because you're multiplying then 2 reals instead of a real and a complex.

The @scalar_rule macro writes rules to use muladd here. @oxinabox does that make a difference?

I have never measured

sethaxen · 2021-10-01T22:17:56Z

src/chainrules.jl

+    return Ω, besselix_pullback
+end
+
+function ChainRulesCore.frule((_, _, _), ::typeof(besseljx), ν::Number, x::Number)


Since everything should be identical for besseljx and besselyx, can you declare them the same with a loop and an @eval?

I could but so far the style in this file is to implement everything explicitly (eg. also derivatives besselj and bessely are implemented separately) and so I sticked to it.

Makes sense.

sethaxen · 2021-10-01T22:24:58Z

src/chainrules.jl

+        ν̄ = ChainRulesCore.@not_implemented(BESSEL_ORDER_INFO)
+        a = (besseljx(ν - 1, x) - besseljx(ν + 1, x)) / 2
+        x̄ = if x isa Real
+            project_x(conj(a) * ΔΩ)


If x is real, then so is a:

Suggested change

project_x(conj(a) * ΔΩ)

project_x(a * ΔΩ)

sethaxen · 2021-10-01T23:24:28Z

src/chainrules.jl

+function ChainRulesCore.frule((_, _, _), ::typeof(besselix), ν::Number, x::Number)
+    Ω = besselix(ν, x)
+    ΔΩ = ChainRulesCore.@not_implemented(BESSEL_ORDER_INFO)
+    return Ω, ΔΩ


My earlier comment might have gotten lost, but since ZeroTangent() * ::NotImplemented is a ZeroTangent, you can implement the frule in such a way that if the AD provides ZeroTangent() for Δν, then the frule does not return a NotImplemented:

Suggested change

function ChainRulesCore.frule((_, _, _), ::typeof(besselix), ν::Number, x::Number)

Ω = besselix(ν, x)

ΔΩ = ChainRulesCore.@not_implemented(BESSEL_ORDER_INFO)

return Ω, ΔΩ

function ChainRulesCore.frule((_, Δν, Δx), ::typeof(besselix), ν::Number, x::Number)

Ω = besselix(ν, x)

∂Ω_∂ν = ChainRulesCore.@not_implemented(BESSEL_ORDER_INFO)

a = (besselix(ν - 1, x) + besselix(ν + 1, x)) / 2

∂Ω = muladd(a, Δx, muladd(-sign(real(x)) * real(Δx), Ω, ∂Ω_∂ν * Δν)

return Ω, ∂Ω

Oh yes, I completely forgot that ZeroTangent() * ::NotImplemented = ZeroTangent(). I'll fix the forward rules!

I noticed that currently it is not possible to test such definitions: JuliaDiff/ChainRulesCore.jl#477

devmotion · 2021-10-02T23:50:25Z

src/chainrules.jl

+    ΔΩ = if ∂Ω_∂ν_Δν isa ChainRulesCore.NotImplemented
+        ∂Ω_∂ν_Δν
+    else


I am not sure if this optimization is useful enough to justify the more verbose and less readable implementation. The optimization is also not performed when the derivative is defined with @scalar_rule.

Yeah I don't think it's necessary. Odds are when this happens it will be because a user actually tried to differentiate wrt the order, and they will instantly realize that's not possible. Ideally the compiler would realize a is unused and not compute it, but that doesn't seem to be the case.

Including it doesn't hurt though. You could always include in a comment the simpler expression so it's easier to follow. We do this in a few places in ChainRules.

devmotion · 2021-10-02T23:51:14Z

src/chainrules.jl

+        a = (besselix(ν - 1, x) + besselix(ν + 1, x)) / 2
+        if Δx isa Real
+            muladd(muladd(-sign(real(x)), Ω, a), Δx, ∂Ω_∂ν_Δν)
+        else
+            muladd(a, Δx, muladd(-sign(real(x)) * Ω, real(Δx), ∂Ω_∂ν_Δν))
+        end


This can't be tested currently and requires JuliaDiff/ChainRulesCore.jl#477 (currently this branch is not reached with neither the default nor the NoTangent() tangent, and finite differencing yields different values (and errors with complex numbers) if a ZeroTangent is specified since then nu is not ignored in finite differencing). With the CRC PR all definitions are tested and tests pass locally.

devmotion · 2021-10-03T23:22:00Z

@sethaxen I think all comments should be addressed. Since JuliaDiff/ChainRulesCore.jl#477 is merged and released also the frules are tested properly now. The tests could be simplified (only one line for frule) if JuliaDiff/ChainRulesTestUtils.jl#222 is addressed; I guess this does not have to block this PR though.

sethaxen

Just a couple minor tweaks to multiply reals before complexes, then this looks ready to merge.

src/chainrules.jl

Co-authored-by: Seth Axen <seth.axen@gmail.com>

sethaxen

LGTM!

devmotion · 2021-10-09T22:02:12Z

Bump 🙂

devmotion · 2021-10-15T22:58:19Z

Some days ago I was added to JuliaMath, so I am able to merge the PR now. I'll wait until beginning of next week, and if nobody objects I will merge the PR then since it was approved by @sethaxen.

devmotion added 2 commits September 29, 2021 01:06

Add derivatives for besselix, besseljx, and besselyx

89e2dc9

Bump version

d1076fc

devmotion mentioned this pull request Oct 1, 2021

Add realdot and imagdot from ChainRules JuliaDiff/ChainRulesCore.jl#474

Closed

sethaxen reviewed Oct 1, 2021

View reviewed changes

Apply suggestions from @sethaxen's review (muladd + optimizations)

c26d35e

sethaxen reviewed Oct 1, 2021

View reviewed changes

devmotion mentioned this pull request Oct 2, 2021

Simplify arithmetic of NotImplemented and treat NoTangent like ZeroTangent JuliaDiff/ChainRulesCore.jl#477

Merged

Improve frules

0af9568

devmotion commented Oct 2, 2021

View reviewed changes

Simplify frules

30a020f

devmotion closed this Oct 3, 2021

devmotion reopened this Oct 3, 2021

sethaxen requested changes Oct 4, 2021

View reviewed changes

src/chainrules.jl Outdated Show resolved Hide resolved

src/chainrules.jl Outdated Show resolved Hide resolved

src/chainrules.jl Outdated Show resolved Hide resolved

Apply suggestions from code review

a9f9f52

Co-authored-by: Seth Axen <seth.axen@gmail.com>

sethaxen approved these changes Oct 4, 2021

View reviewed changes

devmotion merged commit c4b0b83 into JuliaMath:master Oct 18, 2021

devmotion deleted the dw/besselx branch October 18, 2021 14:14

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add derivatives for `besselix`, `besseljx`, and `besselyx` #350

Add derivatives for `besselix`, `besseljx`, and `besselyx` #350

devmotion commented Sep 28, 2021

codecov bot commented Sep 29, 2021 •

edited

Loading

devmotion commented Oct 1, 2021

oxinabox commented Oct 1, 2021

sethaxen commented Oct 1, 2021

devmotion commented Oct 1, 2021

sethaxen commented Oct 1, 2021

devmotion commented Oct 1, 2021

sethaxen left a comment

sethaxen Oct 1, 2021

devmotion Oct 1, 2021

devmotion Oct 1, 2021

sethaxen Oct 1, 2021

sethaxen Oct 1, 2021

sethaxen Oct 1, 2021

oxinabox Oct 2, 2021

sethaxen Oct 1, 2021

devmotion Oct 1, 2021

sethaxen Oct 1, 2021

sethaxen Oct 1, 2021

sethaxen Oct 1, 2021

devmotion Oct 2, 2021

devmotion Oct 2, 2021

devmotion Oct 2, 2021

sethaxen Oct 3, 2021

devmotion Oct 2, 2021 •

edited

Loading

devmotion commented Oct 3, 2021

sethaxen left a comment

sethaxen left a comment

devmotion commented Oct 9, 2021

devmotion commented Oct 15, 2021

Add derivatives for besselix, besseljx, and besselyx #350

Add derivatives for besselix, besseljx, and besselyx #350

Conversation

devmotion commented Sep 28, 2021

codecov bot commented Sep 29, 2021 • edited Loading

Codecov Report

devmotion commented Oct 1, 2021

oxinabox commented Oct 1, 2021

sethaxen commented Oct 1, 2021

devmotion commented Oct 1, 2021

sethaxen commented Oct 1, 2021

devmotion commented Oct 1, 2021

sethaxen left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

devmotion Oct 2, 2021 • edited Loading

Choose a reason for hiding this comment

devmotion commented Oct 3, 2021

sethaxen left a comment

Choose a reason for hiding this comment

sethaxen left a comment

Choose a reason for hiding this comment

devmotion commented Oct 9, 2021

devmotion commented Oct 15, 2021

Add derivatives for `besselix`, `besseljx`, and `besselyx` #350

Add derivatives for `besselix`, `besseljx`, and `besselyx` #350

codecov bot commented Sep 29, 2021 •

edited

Loading

devmotion Oct 2, 2021 •

edited

Loading